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ABSTRACT 

A  horizontal  difference  scheme  that  conserves  both  potential  en- 
strophy  and  energy  for  general  flow  and,  in  addition,  yields  fourth- 
order  accuracy  for  the  advection  of  potential  vorticity  in  case  of  non- 
divergent  flow,  is  derived  for  the  shallow  water  equations  on  the  stag¬ 
gered  grid  as  a  simple  extension  of  the  second-order  potential  enstrophy 
and  energy  conserving  scheme  presented  by  Arakawa  and  Lamb  (1981).  This 
fourth-order  scheme  is  derived  Doth  for  a  Cartesian  grid  and  for  a 
spherical  grid. 

Comparison  by  means  of  numerical  experiments  between  the  newly 
derived  scheme  and  the  second-order  scheme  showed  the  distinct  advantage 
of  the  new  scheme  in  giving  better  development  and  faster  moving  speed  of 
the  law. 
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I.  Introduction 

Recently  Arakawa  and  Lamb  (1981)  (hereafter  AL)  derived  a  second- 
order  potential  enstrophy  and  energy  conserving  scheme  for  the  shallow 
water  equations  in  order  to  improve  the  simulation  of  nonlinear  aspects 
of  the  flow  over  steep  topography.  Their  numerical  experiments  showed 
the  advantages  of  the  scheme  over  the  (potential)  enstrophy  conserving 
scheme  for  horizontal  nondivergent  flow,  not  only  in  suppressing  a 
spurious  energy  cascade  but  also  in  determining  the  overall  flow  regime. 

A  new  horizontal  difference  scheme  presented  in  this  paper  was 
derived  as  a  simple  extension  of  their  scheme  focused  on  the  increase 
of  the  finite  difference  accuracy.  The  scheme,  of  course,  was  designed 
to  conserve  both  potential  enstrophy  and  energy  for  general  flow  and,  in 
addition,  to  give  fourth-order  accuracy  for  the  advection  of  potential 
vorticity  in  case  of  non  divergent  flow;  the  advection  term  leads  to  a 
fourth-order  Jacobian  proposed  by  Arakawa  (1966). 

In  Section  2,  the  sha  low  water  equations  are  presented  and  deri¬ 
vation  method  of  the  fourth-order  scheme  is  outlined.  The  method  is 
much  the  same  as  used  in  AL.  The  derivation  of  the  scheme  is  performed 
in  Section  3.  The  advantage  of  this  fourth-order  scheme  is  demonstrated 
in  Section  4  through  a  comparison,  by  means  of  numerical  time  integra¬ 
tion,  with  the  second-order  scheme  by  AL.  The  appendix  presents  the 
scheme  for  a  spherical  grid  that  can  be  derived  by  analogy  to  the  pro¬ 
cedure  in  Section  3. 

?..  Outline  of  the  derivation  procedure. 

The  governing  differential  equations  for  quasi-static  motion  in  a 
homogeneous  incompressible  fluid  with  a  free  surface  can  be  written  as 

\  V  * 

jy  +  flkxV+  v  (  K+§)  =  o 


(2.1) 


Here  t  is  the  time,  V  the  horizontal  del  operator,  IK  the  vertical 
unit  vector,  V  the  horizontal  velocity,  h  the  vertical  extent  of  a 
fluid  column  above  the  hottom  surface  and  the  (absolute)  potential 


vorticity  q,  the  mass  flux  V*, 

the  kinetic  energy  per  unit  mass 

and  §E  are  defined  by 

s  ( +  ♦  c, )  K." 

(2.3) 

V*2  hV 

(2.4) 

K  a  1/2  V* 

(2.5) 

K+hs) 

(2.6) 

Here  ^  is  the  vorticity,  fc-pxV  ,  }  the  Coriolis  para¬ 
meter,  g  the  gravitational  . deration  and  Ka  the  bottom  surface 

height. 

After  some  multiplications  the  equations  for  the  time  change  of 
total  kinetic  energy  and  potential  energy  in  this  fluid  may  be  ex¬ 
pressed  as: 

^■(KK)  +  V(v'K)  +  '/V5e  -0  (2  7) 


and 


( '/ztf + 9KK5 )+  v  =0 


(2.8) 
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The  summation  of  (2.7)  and  (2.8)  then  yields  a  statement  of  the  con¬ 
servation  of  total  energy. 

It  UCK+fcSHtjMj-O  (2.9) 

where  the  overbar  denotes  the  mean  over  an  infinite  domain  or  a  closed 
domain.  Here  it  should  be  noted  that  the  term  in  (2.1)  involving  q 
does  not  contribute  to  the  change  of  total  kinetic  energy  and  the  last 
term  in  (2.7)  and  (2.8)  cancel  in  giving  (2.9).  These  two  points  must 
be  taken  into  account  in  the  construction  of  the  finite-difference 
scheme. 

The  vorticity  equation  for  this  fluid  motion  may  be  written  in 
the  form 

^-(h?)  +  V(V*?)  =  °  (2.10) 

and  then  we  obtain  the  potential  vorticity  advection  equation. 

+  =  o  (2.11) 

Thus  in  the  absence  of  spatial  gradients  of  q  there  should  be  no  time 
change  of  q.  This  condition  will  be  also  used  to  construct  the  finite- 
difference  scheme. 

Now  hq  times  (2.11)  plus  1/2  q  times  (2.2)  gives  the  equation  for 
time  change  of  potential  enstrophy 

&-(K)*f)  +  r-(Vr*i*)=o 


(2.12) 
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which  leads  to  a  statement  of  the  conservation  of  potential  enstrophy, 

(2.13) 

Our  present  purpose  is  to  derive  such  a  finite-difference  scheme 
for  the  momentum  equation  (2.1)  that  i)  it  is  consistent  with  a 
reasonable  advection  scheme  for  potential  vorticity  advection  equation 
(2.11);  ii)  it  guarantees  conservation  of  total  energy  (2.9) 
and  potential  enstrophy  (2.13)  for  the  general  case  of  divergent  mass 
flux;  and  iii)  it  has  fourth-order  accuracy  for  the  advection  of 
potential  vorticity  when  flow  is  horizontal  and  nondivergent. 

The  derivation  of  the  scheme  is  performed  along  the  following  line. 
First,  we  can  derive  the  same  constraint  as  in  AL  for  the  total 
energy  conservation,  by  the  application  of  requirements  discussed  in 
this  section  to  a  general  difference  scheme  for  (2.1).  Then,  since  the 
scheme  still  retains  a  high  degree  of  freedom,  we  further  impose  the 
following  requirements  to  fix  the  scheme. 

1)  When  q  in  the  finite -difference  analog  of  (2.11)  is  constant  in 
space,  there  is  no  time  change  of  q. 

2)  To  guarantee  conservation  of  potential  enstrophy  the  finite- 
difference  analog  of  (2.13)  must  hold. 

3)  In  case  of  the  nondivergent  mass  flux  advection  scheme  for  (2.1) 
leads  to  the  Arakawa  fourth-order  Jacobian  (Arakawa,  1966). 

4)  The  symmetry  between  the  Cartesian  components  of  the  momentum 
equations  for  the  case  of  a  square  grid  must  be  retained. 


3.  Derivation  of  the  fourth-order  scheme 


The  arrangement  of  the  variables  and  the  indices  on  the  square  grid 
to  be  used  in  this  derivation,  called  the  C  grid,  is  shown  in  Fig.  1. 
Here  u  and  v  are  the  Cartesian  components  of  V  in  x  and  y  directions, 
respectively. 

The  second-order  differencing  for  the  continuity  equation  (2.2) 
can  be  written 


=  0 


(3.1) 


where 


l w t  ,j’vi 

=  L  \tv  ] 


and  we  choose 


I  tv)  _ 

h  i*K,)  “ 


(  h.  )  i*K) 


(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 


here  overbars  -  i  and  -  j  denote  the  arithmetic  average  of  two 
neighboring  points  in  x  and  y  directions  respectively. 

The  general  fourth-order  scheme  for  the  Cartesian  components  of 
the  momentum  equation  (2.1)  may  be  written. 


9 


-  Sip'/iVf*,) +£i*>it}+> i  Unj+h 

+ A;.;+»  td,j+<4  +  oL  (  C  ki-  $ 

-  (k  +  l)z.X/J+!4]  =o  . 

+  X;+t,j+yt  U*+r,j *K  + 

+/s4rA*4j +* 

-(k+$  W,,)-*)  = 


(3.7) 


(3.8) 


where  K  defined  at  the  h  points  is  specified  by 


;i«*J*S4  =  C  >/a  ul  +  1/zll1  j.; 


(3.9) 


Equations  (3.5),  (3.6),  and  (3.9)  are  designed  to  maintain  conservation 
of  total  kinetic  energy  for  the  divergent  mass  flux  (AL).  The  symbols 
a,  8,  y,  6,  e,  <f>,  X  and  y  are  linear  combinations  of  the  q.  X  and  y 
give  additional  generality  to  the  fourth-order  scheme.  Note  that  the 
terms  involving  q  cancel  by  virtue  of  their  form  when  we  derive  the 
equation  for  the  time  change  of  total  kinetic  energy. 

Application  of  (3.7)  and  (3.8)  at  the  points  surrounding  a  ^ 


t 


point  gives  the  finite -difference  vorticity  equation  consistent  with 
this  scheme 


It  ^  V1)’  ““  (f‘-j 

+  Kk,) 

i  V*k,yi  ( fay} i +  )  +'V*'/i'Y  ( ^ 

+  ^ ( ^  f  ^ ) 

-  fcV  )  +  Wi.^W^tE^ 

+  U, (  \  l, )  +  H C,JH%  (Xl,jH  )  ^  ( 3 . 1 C 

where  the  vorticity  change  has  been  expressed  ^  ( klt>£  )  ,•  ■ 


5.  .  =  (  +  +  5  ).-.; 
^ - 

I#  ft 


(3.11 


=  «f  C  U:,j->2-  iti.jA  ~  ^f-Vj  i 


(3.12 


and  h^  is  a  linear  combination  of  h,  as  yet  unspecified. 

Now  we  impose  the  first  requirement,  to  wit,  that  i%/bt 
vanish  when  q  is  formally  set  equal  to  a  constant  on  the  right-hand  side 
of  (3.10),  regardless  of  the  constant.  If  we  write  a,  B,  y,  6,  e ,  <t>,  X, 
and  u  in  general  form  as  linear  combination  of  the  surrounding  q: 


;  v 1  .>  y  » *  A  w*  'A  "•  A 1  VT* 


n 


e)  ,/e\”  JifL  +/g)Y.+  [t\% 


'A 


<l> 


A 


(*) 


ij  “  \/<  ^‘*i+l  **  (/*/  1/C ) 


f3> 


tv> 


(3.13) 


then  when  q  is  formally  set  equal  to  a  constant,  (3.10)  can  be 
written  as 

w<»  i , 

dt  ";'  ~  rf  (MtFKl£^v,-y;v/ip+  CC-E ’X.ulwulyk) 

+  (B-F)(U’v,-^-)4  CP-E)(  it* -,,)*) 

^  (C~  F  )(  -  IT t-x,]-i)  +  ^  ^  E )  { lAt,j-'/t-  LLi-yj-K  ) 

-KP+F K U?*.; -  V*fc>;W )  +  ( &  +  F )  (  U»i,j-k - 

+  M  ( -  Vw-.; > 

-  L  (u-yy^u^-u^-.y.-ic' ,„*)] 


(3.14) 


(  K  ) 

where  ^  s£c<  ^ 


oo 


,  etc. 


* 
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In  the  case  of  a  square  grid,  for  simplicity  and  geometrical 
symmetry,  is  defined  as 

C4>  I  i 

=  4  (  +  ) 

Then  from  (3.1)  and  (3.15)  we  can  write 


(3.15) 


hi.j  -  -  *^j-  £  (  ~  } 

■*“  ( ,;v —  + 

*  (  ^t-4,  j  -  ^  iLi.'j-K  ~~  O 

+  C  ~  +  ILw.y/i  ~  W**'.  j-h.  ^  3 

Comparison  of  (3.14)  and  (3.16)  yields  the  constraints, 
A=B=C=D3  1/4,  E=F  =  L  =  M  =  0 


(3.16) 


(3.17) 


For  the  requirements  2),  that  potential  enstrophy  be  conserved 
in  finite  difference  form,  we  formally  rewrite  the  vorticity  equation  as: 


+  +  b*';fo=0 


(3.18) 


*  ★ 


where  d/.j  and  b.\j  are  linear  combination  of  U  and  V  . 

Then  according  to  AL  the  necessary  and  sufficient  conditions  for  the 
conservation  of  potential  enstrophy  in  divergent  flow  can  be  formulated 


&L,j  *  t+k/  J+j  “  > )  *j':  *  ’ j 

b<*.j  ~  ^ 


(3.19) 


(3.20) 
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To  impose  these  constraints  we  use  (3.13)  to  rewrite  (3.10) 
explicitly  in  the  form  giyen  by  (3J8)  and  obta'n  the  following? 


b,j  =  a-  f)  *  (f  f)  ^ . ♦  ( 6“v  r>  < 


,,»)  Of  .* 

lA  MX.;  ~J* 


-  (f+£“4l4-(/-  f)  U.*;.!4  4  «V') 


(3.21) 


fcy  =  (o("-SV’>  +  (  tf-Kj  > 

+  (£  J  $%)-)  +  )f  ~  W*4  A 

-  f'lrV-  ( (SVU'K* 

-((AeV*.)-*-  ltVx. 

=  U°-A/Ai4;  ♦ 

-  /a'*'  /i°'i  ,k*  _  /  r*'  *.<*'  0«, .  .*  /,'!■;  ,<li  ,  * 

Y  “  t  C  A  —  <T  -t-X  )lCi,ji!i  +  (S  -t  >U.^,j.i/x 

'  oJ  ”  |3  +-A.  )  t  (<X  -  C  ) 


Ai.jn.jn  =  (=i'-t;/")^4((3  n'-ry  ■>***.; 

~  ( (A  f>T.v-  ( A  Ou'^-CS-irV^r.;.,. 

'  (6'"-  t'">  A,;,V, 


+ < £  n+  f  ^  4  C  J  -  <f “) U;-x,;-i  - j., ,j, i 

/r  .  .i.  ^  i/.  .a  . 


jr>  l|>  H;  . 


U,T-  (AX')  U*n  + 


(3.22) 


a;, =  -  rvk-r 

-/v^vj  -  ( r-  n  k,v;* 

-  ( r  *"')  u\r%  -  fua^h  +xv,.,.* 

Ot;: Hr*  -  «lX*  t  (3'^v;  +  ( f)  A,j-i 
+  5 mvi,.r, - t  «*%  n  ici.j-'y, 

~~  ( d*- (3<V>)  +  fP'lLu,)* \  i  A. 

Cil.j  ■  H}»  Sr  ~  £  if  ^  “  ft  “dv  V'Zb'p, 

-  ( *V,  -y"»U:}  -  < ^ 

+(f-nuLj+tf/uLM  +j?ui}* 

fcj :  !*'rj-i  -  ci  +  6U^Vi.j  +  c  £*+  <f>  h*w,y/ 

^  r*v< 

~  +  £  )  M'i*i,j-'/l  “  W  -p  )  Ustya  +K  Ustyik 

d,i.y  in ,~)  ~  ~  J*'  ^;*<\}  f  U* •■«,/*'*  f  U-^y/x 

r'?'N*  J9,\r*  /i‘7\r*  i*\c* 

O  tXVt;,v. 

fl;.j:i.i-i  =  oi^iTi'n.j  *  p^V**,)  ♦  S<Vf**.;-i  + 

-U','-(3"’)U.!;-!i  +X'l£j-,* 

.  <u  ^  .(vi  if 

dij  :  (ti,)4|  -  -  0  ,  U;,yin,}-I  *  <* 

CU)  •  i+»,j+2  =  “  If  }li'<+\)*'/x  f  (X;,)'-i-*,}-l  —  o(  Ut-Kr'/i 

u,  *  c^.i* 

CLi.;:^l'rL=  “  P  U'W.'rh  >  0  U,<-i,}+h 

Ar.j  :  1‘^jtZ  =  £  »  &y*  I*1’)-2-*  ~  P  Ul+>,j- 


Applying  first  the  constraint  (3.20)  fpr  arbitrary  U  and  V*,  and 


simplifying, using  (3.17),  requires 


sw=  r=a 

r+<r=  r-fw* 

,T)  x(l)  if/ 

o(  +  <P  =0  -</>  =!/? 

r-?‘"=  r'- £'*'=(/? 

e"-  oT’-tg^.l/i 

X”  -  a'3'  =  o 


(3.23) 


To  apply  the  constraint  (3.19),  each  CLit  j  i-£'  j- >' 
required  to  equal  -  ft.-;  :  )■»>'  when  all  pairs  of 

indices  appearing  in  the  former  are  incremented  by  (£'yj')  .  This 
gives  the  following  conditions: 


(3.24) 


or  =  a**.  «-“= r=  r  -  5=  * 

A*"--/" 

oT’«  <*‘w  =  -  -  xT 

cl*'  -  (3“'  -  f*’-  «*'»  Vi*  -A“- yt/'’ 


For  the  requirement  (3),  to  wit,  the  reduction  of  the  non-divergent 
mass-flux  advection  to  the  fourth-order  Jacobian,  we  rewrite  (3.22)  for 
the  case  of  no  mass  flux  divergent  introducing  a  stream  function.  By 
equating  the  resultant  with  the  fourth-order  approximation  of  the 
Jacobian  (Arakawa,  1966),  we  can  determine  the  following  coefficients 
(See  Appendix  A) 


3.25 


and  thus,  from  (3.24) 


tt>  <*> 

0/  =  (3  ’ 

cr3’=(r 
dh)-  f 


X“'=SM- 

/” = r  ■ 


'/e 

-  l/2* 
-  ‘/^f 


Now  (3.23)  and  (3.24)  can  be  combined  to  express  unspecified 
coefficients  in  terms  of  fc"',  £**’,  f  ”  and  f(,>  : 

a  "  =  Vt-l",  oL^/9-r',  d“'=-l 

e'”=  >/i++£!,i/ eT- •//-«: -4’"’ 

S"  V?  -t  £'3’-  f"  b*'~ .  /*+  +  ^ 

£(I  =  1/24 £“"^ ¥*- ~ Viz'^"’j 

Finally  the  requirement  of  symmetry  between  the  Cartesian 
components  of  the  momentum  equations  for  the  case  of  a  square  grid  yields 


(3.28) 


(3.26) 


(3.27) 


.(f) 

(z;_ 

.<r; 

(i> 

1 ^ 

p 
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If  we  now  choose,  simply  for  reasons  of  symmetry  between  6  and  $ 
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the  terms  in  (3,13)  are  completely  specified  as; 
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When  we  disregard  additional  coefficients,  X  and  p  in  the  momentum 
equations,  it  is  easily  shown  from  (3.23)  -  (3.30)  that  the  final  forms 
of  a,  3,  y,  6,  e  and  <j>  coincide  with  those  of  the  second-order  scheme 
presented  by  Arakawa  and  Lamb  (1981). 

In  summary,  use  of  (3.5),  (3.6),  (3.9),  (3.11),  (3.12),  (3.15) 
and  (3.31)  in  (3.1)  -  (3.4),  (3.7)  and  (3.8)  give  a  fourth-order 
potential  enstrophy  and  energy  conserving  scheme  for  the  shallow  water 
equations. 


4.  Numerical  test  of  the  fourth-order  scheme 

In  this  section  the  potential  superiority  of  the  newly  derived  scheme 
is  examined  by  means  of  the  numerical  experiments  comparing  its  results 
with  those  of  the  second-order  scheme  of  Arakawa  and  Lamb.  For  this 
purpose,  a  two-layer,  homogeneous,  incompressible  fluid  system  model  was 
used.  Quasi-static  motion  in  the  each  layer  of  this  fluid  system  may  be 
described  by  (2.1)  -  (2.6),  but  replacing  h  in  (2.6)  by  rka  +  lu; 
here  (and  hereafter)  subscripts  1  and  2  denote  the  lower  and  upper  layers, 
respectively;  and  r,  equal  to  1  in  the  upper  layer  and  to  the  ratio  of 
density  of  two  layers  in  the  lower  layer,  expresses  the  effect 

of  stratification.  In  this  fluid  system,  the  summation  of  total  energy 
for  two  layers  is  conserved,  and  mass  and  potential  enstrophy  are  conserved 
in  the  each  layer. 

In  the  numerical  experiments,  X  and  Y  are  directed  toward  the  east  and 
the  north,  respectively.  The  computational  domain  is  confined  in  a  channel 
bounded  by  two  parallel  rigid  wall  at  y  =  0  and  y  =  W  =  7500  km;  and  by 
X  =  0  and  X  =  L  =  13500  km,  where  cyclic  boundary  conditions  are  applied. 

The  mean  heights  of  each  layer  top,  H.j  and  H^,  are  chosen  to  be  5  km  and 
10  km,  respectively  and  grid  size  d  =  500  km.  The  bottom  topography  is  a 
circular  mountain,  placed  at  the  center  of  the  domain,  of  circular  and 
parabolic  cross  section,  and  expressed  by 

hs  *  hm  (1-D2/R2)  for  0  4  D  <  R 

and  h$  =  0  for  R  <  D 

here  0  is  the  distance  measured  from  the  center  of  the  mountain,  R  =  3  d  so 
that  the  bottom  diameter  of  the  mountain  is  3000  km  and  the  maximum  height, 
h  ,  of  3  km  is  chosen  (see  Figure  2).  The  ratio  of  density  r  =  2/3,  and 
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_2 

the  acceleration  of  gravity  g  =  9.8  m  s  ,  are  selected.  Experiments 

were  performed  with  two  different  Coriolis  parameters:  one  is  a 

-4  -1 

constant  f  =  f  =  10  s  (f-plane)  and  the  other  is  variable 
f  =  fQ  +  8  (y-W/2)  with  8  *=  1.62  x  10'^nf^  s"^  (8-plane) .  At  the  north 
and  south  boundaries,  the  rigid  wall  conditions,  v  =  0,  and  a  computa¬ 
tional  boundary  condition,^  =  0,  are  applied.  The  variable  outside 
the  boundaries  involved  in  the  fourth-order  scheme  are  expressed  by 
those  of  interior  points  using  the  conditions  that  u,  h,  f  are  sym¬ 
metric,  and  v  is  anti -symmetric,  with  respect  to  the  boundaries.  The 
preliminary  numerical  integration  with  these  conditions  showed  that 
they  worked  well  within  the  acceptable  accuracies  for  conservative 
quantities:  0.04%  decrease  in  mass,  0.08%  decrease  in  total  energy, 
0.04%  and  0.03%  increases  in  potential  enstrophy  for  the  lower  and 
upper  layers,  respectively,  after  10  days.  However,  a  certain  propor¬ 
tion  of  the  accuracy  may  be  attributable  to  round-off  error  from  the 
use  of  single  precision  in  the  code.  As  initial  wind  fields,  uniform 
westerly  flows  of  20  ms-^  and  30  ms"^  are  prescribed  for  the  lower  and 
upper  layers,  respectively.  The  height  fields  are  then  obtained  from 
the  steady  state  solutions  of  the  equations  for  the  two  layer  fluid 
system  in  the  absence  of  a  mountain.  The  Matsuno  scheme  is  used  first 
for  the  initial  time  step,  and  then  once  in  every  five  leapfrog  time 
steps.  A  time  interval  of  eight  minutes  is  applied  for  all  experi¬ 
ments,  though  the  second-order  scheme  may  allow  use  of  a  longer  time 
interval  than  that  acceptable  for  the  fourth-order  scheme  (Gerrity 
et  al . ,  1972).  This  time  integration  scheme  is  the  same  as  that  used 
in  the  UCLA-GCM,  for  reasons  of  convenience.  The  model  time  of 
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integration  was  20  days. 

Fig.  3  shows  the  relative  time  changes  in  potential  enstrophy,  total 
energy  and  mass  to  their  initial  values. 

For  all  cases,  total  energy  and  total  mass  monotonously  decrease 
almost  at  the  same  rate  up  to  20  days  (0.15%  and  0.08%,  respectively). 

These  rates  of  decrease  are  much  smaller  than  those  resulting  from  the 
fourth-order  scheme  presented  by  Navan  and  Riphagen  (1979)  (1.44%  ^r  0.99% 
for  total  energy  and  0.54%  or  0.42%  for  mass,  after  4  days).  For  potential 
enstrophy,  the  fourth-order  scheme  gives  better  conservation  than  the  second- 
order  scheme  does. 

Figs.  4-23  show  the  24-hour  changes  of  the  wind  vector  based  on  u1  and 
v1  at  h  points  and  the  differential  height  of  the  lower  layer  (height  of 
interface)  h-j  -H^ ,  for  a  period  up  to  20  days,  for  the  fourth-  and  second- 
order  scheme  on  the  f-plane.  Comparison  of  these  results  shows  that  the 
fourth-order  scheme  yields  a  faster  translating  speed  of  the  low  (inter¬ 
facial  depression)  produced  by  the  mountain  (the  effect  of  truncation  error 
in  general  being  to  reduce  propagation  speeds).  The  five-day  average  speed 
of  the  lows  in  the  fourtn-  and  second-order  schemes  are  19.1  ms~^  and 
16.8  ms-1,  respectively,  so  that  the  difference  of  distance  of  low  centers 
reaches  about  1000  km  after  6  days.  The  passing  event  over  the  mountain  at 
Day  8  and  the  subsequent  reformation  of  the  traveling  low  on  the  lee  side, 
consistent  with  potential  vorticity  conservation  are  dramatically  reproduced 
in  the  fourth-order  scheme,  but  not  so  distinctly  done  in  the  second-order 
scheme.  In  the  second  passing  event  after  Day  15,  which  exhibited  somewhat 
different  behavior  from  the  earlier  one,  the  difference  of  the  results  of  two 
schemes  is  also  clearly  seen.  The  deep  traveling  lows  in  the  fourth-order 
scheme  pass  just  north  of  the  mountain,  preserving  their  forms  and  'n- 
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tensities,  and  also  producing  the  jet  stream  near  the  boundary  associated 
with  the  stationary  high  over  the  mountain(Day  16),  while  in  the  second- 
order  scheme  the  weak  low  is  absorbed  by  the  relatively  strong  high  as  it 
approaches. 

The  results  of  the  6-plane  computations  are  shown  in  Figs.  24-43. 

In  this  case  the  initial  height  fields  are  not  of  course  identical  to  those 
in  the  f-plane,  but  the  differences  are  so  small  as  to  be  negligible.  Thus 
we  can  determine  the  effect  of  6,  comparing  the  results  in  this  case  with 
those  in  the  f-plane  case.  The  evolution  of  whole  pattern  in  the  6-plane  is 
much  more  complex  and  the  intensities  of  the  high  and  the  low  are  greater 
than  those  in  the  f-plane  case.  One  of  the  remarkable  features  is  the 
occurrence  of  the  splitting  of  the  lee  side  low  into  a  stationary  low  and  a 
low  traveling  toward  the  northeast  of  the  mountain,  never  seen  in  the  f-plane 
case.  Note  that  the  occurrence  time  of  the  splitting  in  the  fourth-order 
scheme  is  one  day  earlier  than  that  in  the  second-order  scheme  (see  Day  4 
and  Day  5) .  Here  we  can  see  the  same  advantages  of  the  fourth-order  scheme 
as  noted  in  the  previous  f-plane  computations.  After  6  days,  the  difference 
of  distance  of  the  traveling  low  centers  has  increased  to  2000  km,  twice  as 
much  as  that  in  the  f-plane  case.  Later,  the  difference  gradually  decreases 
because  the  low  in  the  second-order  scheme  is  catching  up,  while  the  low  in 
the  fourth-order  scheme  is  blocked  by  the  strong  blocking  high  over  the 
mountain.  The  traveling  lows  can  not  easily  pass  over  the  mountain,  but 
go  around  it  until  they  combine  with  another  low  in  the  south,  associated 
with  the  planetary  wave.  This  splitting  and  the  associated  blocking 
phenomenon  are  not  produced  in  the  similar  experiment  by  Kasahara  (1966). 

The  motion  in  the  well-developed  strong  high  over  the  mountain  occasionally 
resembles  a  Taylor  column,  with  a  stagnation  area  just  south  of  the  mountain, 
although  this  effect  is  confined  to  the  lower  layer. 


22 


5.  Summary  and  further  comments 

A  horizontal  difference  scheme  has  been  derived  for  the  shallow 
water  equations,  conserving  both  potential  enstrophy  and  total  energy 
for  general  flow,  and,  in  addition,  yielding  the  fourth-order 
accuracy  in  the  potential  vorticity  advection  (in  case  of  horizontal 
nondivergent  flow),  in  the  presence  of  bottom  topography.  This  de¬ 
velopment  may  be  considered  as  a  straightforward  (if  complicated) 
extension  of  the  second-order  potential -enstrophy  and  energy  con¬ 
serving  scheme  of  Arakawa  and  Lamb  (1981).  Comparisons  by  means 
of  numerical  experiments  using  the  two-layer  model  with  the  second 
order  scheme  by  Arakawa  and  Lamb  showed  the  advantages  of  t  v  »»<wly 
derived  scheme  in  better  development,  faster  traveling  speeds  of  the 
lows  produced  by  the  mountain,  and  better  conservation  of  potential 
enstrophy.  The  increase  of  computation  time  by  the  fourth-order 
scheme  is  less  than  10%.  On  the  present  model,  the  calculation  code 
of  the  advection  terms  occupy  the  major  part  of  the  whole  code,  so 
that  the  increase  in  computation  time  may  become  negligibly  small 
in  more  complex  GCM  or  NWP  models.  The  new  scheme  has  been  already 
incorpolated  into  the  current  UCLA-GCM  with  the  SICK-proof  kinetic 
energy  for  the  scheme  defined  by 

- 1 

Ki*  =  f  '/a  U)+yz  +  Vza  (  ^ j+  Vz + )z  3  it  y* 

+  (  '/3  u\ty2  +  '/24  c  lTc+j/i+  (5.1) 

This  has  practically  eliminated  the  linear  computational  instability 
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Of  meridionally  propagating  inertia-gravity  waves  (Arakawa  and  Lamb 
(1981)),  and  has  given  good  simulation  results,  GISS  also  has  used 
the  scheme  for  their  climate  model. 
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APPENDIX  A 

The  fourth’-order  scheme  for  the  case  of  no  mass  flux  divergence 

The  finite-difference  Arakawa  fourth-order  Jacobian  J^4^  (1966), 
which  originally  represents  the  advection  term  in  the  vorticity  equation 
and  guarantees  the  enstrophy  and  kinetic  energy  conservation  for  hori¬ 
zontal  nondivergent  flow,  can  be  directly  applied  to  the  present  case 
of  potential  vorticity  and  potential  enstrophy  because  of  its  universal 
form.  Thus  the  schemes  for  the  momentum  equations  (3.7)  and  (3.8)  con¬ 
serve  potential  enstrophy  and  have  fourth-order  accuracy,  if  the  finite- 
difference  potential  vorticity  equation  (3.10)  derived  from  it  reduces 
to  this  Jacobian  form. 

The  Arakawa  fourth-order  Jacobian  for  the  present  case  can  be 
written  as: 


L 

-  (  t  .  )  ( t i  ) 

+  (fkj  +fV"  fUr'&rrXfo* 

+  < ( %r  Kl, ;)($.■.;<  Ur) 

*  ( - < %-%-,)%  +u/r, )J 

U-ri)  '(%’.]■> r  i:*;-) Ul-i) 

-  (  -  ft,. j  )  %  Hi-  +  C  fw,  ■  -  f*rl)  % 

-  c  + m, 3 


(A.l) 
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where  •/'*  denotes  a  stream  function  for  the  mass  flux. 
With  the  definition  of  9*  as 


d'  ( f*  -  %») 


the  right-hand  side  of  (3.10)  can  be  rewritten  as 


(A. 2) 
(A. 3) 


(-otq-4 +S,,j*x 

+  4?i,j  ( f ; ,  £;-x,  •*+/'  '0  +■/*  '-';P 

+  X-'i+rhu.  i‘H~  1 

+  •/(«,;-»  t-  &'%•/-£  ^ 

+  ‘fiT;-»  C-Atf-)+f4JXM«> 

+  'fA.J  ^  ~ Ml~[']  1  Vi«,;  (/<•«. j)  ^  (A-4) 


Equation  (A. 4)  with  the  right-hand  side  of  (A.l)  gives  the  forms  of  X 
and  y  and  four  independent  constraints  on  a,  6,  y,  6  ,  c  and  $  as 


functions  of  q: 
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After  some  algebraic  manipulations  similar  to  those  in  AL, 
we  obtain  the  final  forms  as; 


f  16  0&  &•)+' +  ^ ?Z'j + &''}+ ' } 

+P*h)+0) 

~  z(Z+ }+'/i 

^ —  ~  ^1%,  J+%  ^  ~fi  \J^  il+ft.j  )  ~3  ^ /tJ  i+hj+l  ) 

+  %i-vr i  (fa*.]+i +  ?) 

“zc*+f  )*•*,/<& 

^  =  Ci-4,  j*  +  ,7  fa  ti,j  +* ( faj+k/t)  -3  +£'«,;  ? 

+  <HjW  +?HJta-  (ft-J  J  +fi-t/}-r+frfn + 

+*(*+♦)*%.!* 

+  ?£+.,Jt2+ftt2j+r  (fa,]H+tt'j*2+  j)) 

Xt,J+l=  5*  &-b>+l  ) 


«•,)  (  i W,/-l  fit4*,)*1) 


(A. 6) 


where  the  new  variable  C  is  defined  by 


r  -  I 


(A. 7) 


APPENDIX  B 


The  fourth’-order  scheme  for  the  shallow  water  equations  on  a 
spherical  grid. 

1.  The  governing  equations  in  spherical  coordinates. 

The  shallow  water  equations  (2.1)  and  (2.2)  on  a  spherical  grid 
may  be  described  as: 

h(w)-^h(K^)=0  (E 

ft(£)+8^+!^K+f)  =  o  (b 

h(h)+h,^r^k^)=0  (B 

where  m  and  n  are  metric  factors  in  the  spherical  coordinates  %  =*  X 
(longitude)  and  ^  (latitude)  expressed  by  1/m  =  a  cose/* 

and  1/n  =  a  (the  radius  of  the  earth).  U  and  V  are  the  component  of  V 
in  \  and  ,  respectively  and  q,  K  and  $  are  defined  in  Section  2. 
The  vorticity  ^  —  IK  V'XV  can  expressed  as 

«nr>  .£_>u  i  (B 

Then  the  equation  for  the  time  change  of  kinetic  energy  can 
be  obtained  as 
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2.  The  finite-difference  scheme  for  interior  points  of  the  spherical 
grid. 

A  portion  of  the  spherical  grid  with  the  variables  staggered  as 
in  the  C  grid  and  the  indices  (t.j)  centered  at  a  q  point  is  shown  in 
Fig.  B1 .  Here  4^  and  4^  are  constant  grid  intervals  in  ^ 
and  \  ,  respectively,  and  m  and  n  are  assumed  to  vary  only  in  j. 

For  the  continuity  equation  (B.3)  multiplied  by  4^4>£  »  the 

following  finite  difference  form  is  chosen: 

(B.6) 

(B.7) 

(B.8) 

Here  4^4*2  is  area  of  the  stippled  region 

in  Fig.  B1  and  h(u),  h(v)  are  as-yet-unspecified  functions  of  h. 

Ignoring  pressure  gradient  forces,  the  remaining  terms  in  (B.l) 
and  (B.2)  can  be  represented  as  in  the  square  grid  case. 
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but  U  and  V  are  now  defined  by  (B.8).  The  coefficients  a,  8,  y,  <S,  e, 

( p ,  X  and  y  are  again  functions  of  q  to  be  determined,  where  fl..a  {*? 

and  h^  is  an  as-yet-unspecified  function  of  h.  The  finite-difference 
form  for  ^  chosen  is 


Conn 


(B.ll) 


As  in  the  square  grid  case,  h^  and  h^  are  chosen  as 


(B. 12) 


Then  it  is  easily  shown  from  the  finite-difference  analog  of  (B.5)  that 


we  must  specify 


IW* = (~^  [|  u‘‘  + 


(B. 13) 


to  maintain  conservation  of  total  kinetic  energy  for  divergent  mass  flux. 


Since  tne  procedure  to  determine  a,,  y»  <5,  e,  <f>,  \  and  y  as 
functions  of  q  Is  identical  to  that  presented  for  the  square  grid,  we 
can  directly  use  the  results  of  (3.31)  with  consistently  determined  h 
in  q  for  the  spherical  grid  (AL ) ; 
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3.  Modification  of  the  scheme  at  the  pole 

The  arrangement  of  the  variables  and  the  indices  near  the  North  Pole 
(NP)  is  shown  in  Fig.  B2. 

The  continuity  equation  at  j  =  NP-1/2  in  the  grid  system  can  be 
expressed  as 

H ws.NP-t*4  U;*  i  vi  - 16*  .  w-'/z  -  lrt+y2 ,np-i  =0  ( b.  1 6 ) 

The  momentum  equations  for  u  at  J  =  NP-1/2  and  V  at  j  =  NP-1 
can  be  written 


P-'/z  $-\NH 

+  ii+VisW-t/i  lli*i,NP-t/2'‘Zt-yl,Nf~Vt  U-Hf/r-yz  H 

+  [  3  =  o 

h*<, Nf~ >/i  + (8 l£ )i,NF-yz 4  W v? )i.Np - V2 
4  (  S  16*  )i*\NF-i/i  ~  fi+VijNr-l  l)Z*Vi,*P-i  6^tH,NM 

~  JAi,nr- 1  irZyirHP-l  4  t  K iWx,NF-Vi.~  Ki*Vi,NP-Vj.  3  =  0 


(B.17) 


(B.18) 


where  the  primed  variables  are  to  be  determined  and  the  others  are 
defined  as  in  (3.31). 


J 


At  the  North  Pole,  we  define  q  through  the  circulation  theorem 


where 


(B.19) 
(B.20) 
(B.21 ) 


The  factor  in  (8.21)  represents  the  area 

of  the  hatched  region  in  Fig.  B2,  and  IMAX  is  the  total  number  of  grid 
points  in  the  ^  direction. 

For  the  conservation  of  total  kinetic  energy  in  the  divergent 
mass  flux,  k‘  at  j  ■=  NP-1/2  must  be  given  by 


'  -  1  (WOmm 

1  A%A>i 


(win  )nr-vt 


Cm ) 


(B.22) 


After  the  application  of  requirement  on  potential  vorticity 
advection  and  potential  enstrophy  conservation  to  (B.16)-(B.21),  we 
obtain  the  constraint  on  the  area-weighting  factors,  the  form  of  hp^q 
and  y'.  <$'  and  c'  as  functions  of  q. 
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Cl nh)itp-i 


Crt-fl  )vp-  Vz  C  win  )wp-  v» 


(B.23) 


rt)  j  IM4X 

"P  ”  JM/JX  k*Mr-*k 

fljL  _  i  _ASii 

(wftV  2  Orrw.)Nr.yl 


(B.24) 


(B.25) 


^ i,NP~'/i  =  L  3  +  2  Np./  +  3  $  I(  |  ^  np-  (  ) 

5c,NP~Yz  =  [3^+3  ?l,AfP-l  +  2  t*t4,,/vP-/ 

£  i+/4,NP-/j  =  [  2  $\p  ~~  1  —  S  i»'/NP-»  3 

^i,wr-}i  —  ^4  [3  £;,vP-(  +  ' $*p] 

/^•'♦'/NP-I  *  l!4  ^  %r-*',/»P-Z  “  ftvp} 


(B.26) 


(B.27) 
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The  proceudre  used  for  the  North  Pole  can  be  directly  applied  to 
the  South  Pole  (SP)  taking  Into  account  geometrical  symmetry,  so  that 
we  give  only  the  results  here. 

.  The  continuity  equation  at  j  =  SP+1/2 


^  If  y  y 

M  ,  _  AUt  . 

Ht^SrA  "  0*irtW,  hl*>"** 


(B.28) 


The  momentum  equations  for  u  at  J  =  SP+1/2  and  v  at  j  =  SP  +  1 


dfrv  °^/*P**i  ^Wi,spM 

/  %  /  $  * 

+  ILitSf+’A 

"*  C  Ki-\sr*yt^  ~  ® 

bi  ^YL~  X^,sp+/5+  t  ^  h'WK  t  CSli  )i,sr*>A  ^  U  u*  hsMi 

~  yUi,Sf*(V}.yllsr*>  *  (K.H*,sr*K  “  K  3  *  O 


(B.29) 


(B.30) 


( B. 31 ) 
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and 


/ 

«u 


» — 


~i£  U  i»i,w + 3  i+3  ftp  -  ft,£r*2  -  ft-,/SP, ,  J 

P«,Sp*Ji-  24  (3  ft-ySF+t  +<2  ft-l/Sf4l'<'3ftp“^»/SfHI'"?t',5f+2  ) 
-24  f  +  ft,sp+l  “  2  ftp  3 


^ ;  ,sp+  %.  = 
Si,  Sp+3^  = 
/^?+vsp+i  = 


J 


( B .  32 ) 


1 

if  t  ^  &.3F+2  +  +  2  ftw,3p«  •+  3  ?t,sp+'  ”?*fp"?znsp+f3 

24  [  f<>'/5P*2  +  2ft/5ft2  +3 
*4  t  J 


Potential  vorticity  and  vorticity  are  defined  as 


2  =  hrj-5*t 

bs?  Ltt* 

Asp 

(B. 34) 

e  —  X^7*4*\ 
bsp  4Jr^J (  ^ 

(B.35) 

=  IMAX^m 
r  flWSp 

(B. 36) 

(f)  |  IMAX 

K,f  IMAX 

(B.37) 

*W  l  AiAt, 

C'Wl7t)Sp  2 

(B. 38) 

1  (•  4*44  1 

t/7in.)spf,  2  Lflrin^p^  OunJspp-H.  J 

(B. 39) 
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LEGEND 

Fig.  1  The  staggering  of  the  variable  based  on  the  C  grid. 

Fig.  2  A  horizontal  domain  with  a  circular  mountain  and  a 
vertical -section  view  of  the  two-layer  fluid. 

Fig.  3  Time  change  of  the  conservative  quantities  to  their 

initial  values  expressed  by  the  percentage. f  :  constant 
f  case,  3p:  variable  f  case,  the  solid  line:  the 
fourth-order  scheme,  the  dashed  line:  the  second-order 
scheme. 

Fig. 4-23  Results  of  the  time  integration  for  the  constant  f  case. 

The  arrows  and  contours  represent  the  wind  vectors  and 
the  differential  heights  of  the  lower  layer,  respectively 

Fig. 24-43  As  in  Fig.  4-23,  except  for  the  variable  f  case  (e-plane) 

Fig .81  A  portion  of  the  spt,  rical  grid.  The  area  of  the 

dotted  region  is  represented  by  (A^d^/mn.  v*  . 

Fig.B2  The  spherical  grid  near  the  North  Pole.  The  arees  of 
the  dotted  and  hatched  regions  respectively  represent 
)MP  and  (A^/mn  Juf.yj 
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Second-  and  Fourth-order  Schemes 
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